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Abstract 


Forward-backward sweep method (FBSM) is an indirect numerical method 
used for solving optimal control problems, in which the differential equation 
arising from this method is solved by the Pontryagin’s maximum principle. 
In this paper, a set of hybrid methods based on explicit 6th-order Runge— 
Kutta method is presented for the FBSM solution of optimal control prob- 
lems. Order of truncation error, stability region, and numerical results of 
the new hybrid methods were compared with those of the 6th-order Runge— 
Kutta method. Numerical results show that new hybrid methods are more 
accurate than the 6th-order Runge-Kutta method and that their stability 
regions are also wider than that of the 6th-order Runge-Kutta method. 
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1 Introduction 


Numerical methods used for solving optimal control problems (OCPs) are 
generally divided into two categories, direct and indirect methods. Indirect 
methods solve an OCP numerically based on the Pontryagin’s maximum 
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principle. The Pontryagin’s maximum principle was proposed in 1954 by 
Pontryagin, a Russian mathematician. In the indirect method, an OCP is 
converted into a two-point boundary-value problem (TPBVP). The forward- 
backward sweep method (FBSM) is one of indirect numerical methods, which 
was first proposed in 2007 in a book entitled as “Optimal Control Applied 
to Biological Models” written by Lenhart and Workman [13]. In 2009, Sil- 
veira et al. [23] suggested six methods for classifying skin lesions in medical 
images and concluded that the FBSM performs better than the other meth- 
ods. In 2012, McAsey, Moua, and Han [16] proved the convergence of the 
FBSM. In 2015, Moualeu et al. [17] used the FBSM to treat and control a 
type of tuberculosis with unknown cases in Comeroon. In 2015, Rose in his 
thesis [20] reported that the FBSM is more accurate than the direct shooting 
method and optimization method of MATLAB . In 2015, Sana et al. [22] used 
trapezoidal and Euler methods, instead of Runge-Kutta method through the 
FBSM to solve an OCP, compared them with the 4th-order Runge-Kutta 
method and concluded that trapezoidal and Runge-Kutta methods have the 
same performance in solving the OCP. Indeed, the performance of these 
two methods was improved by increasing step length compared to the Euler 
method; see [21]. In 2017, Lhous et al. [15] proposed a discrete mathematical 
model and optimal control to reduce the divorce rate. They used the Pon- 
tryagin’s maximum principle and FBSM. They informed the people of the 
community about advantages of marriage and disadvantages of divorce, and 
in this way, they were able to reduce the divorce rate. In 2018, Kheiri and 
Jafari [11] proposed a general formulation for a fractional optimal control 
problem (FOCP), in which state and co-state equations are given in terms 
of left fractional derivatives. They used an improved FBSM, by the Adams- 
type method to solve the FOCP . In 2019, Duran, Candelo, and Ortiz [3] 
used a modified FBSM for reconfiguring unbalanced distribution systems. In 
2019, Kongjeen et al. [12] proposed a modified FBSM for analyzing electrical 
charge of microgrids. In 2020, Ameen, Hidan, and Mostefaoui [1] proposed a 
mathematical model for studying the relationship between fish consumption 
and the prevalence of chronic heart disease (CHD). They used an improved 
FBSM based on a predictor-corrector method and concluded that eating fish 
reduces the risk of CHD and its mortality. In 2020, Bhih et al. [2] pro- 
posed a new model of rumor on social media and determined three optimal 
controls theoretically minimizing the number of spreader users, fake pages, 
and related costs. They used the FBSM to solve their optimization system 
in a duplicate process. In 2021, Ebadi et al. [9, 8] presented hybrid meth- 
ods to numerically solve the OCP by the FBSM. They concluded that their 
proposed methods are more accurate than the FBSM in the presence of the 
explicit Runge-Kutta methods. 


In this work, we present new hybrid methods for the FBSM solution of 
OCP. The paper is organized as follows: In Section 2, new methods and their 
order of truncation errors are presented. Stability analysis of the methods 
discussed in Section 3. We review the OCP, FBSM, and its convergence in 
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Sections 4— 6. The results and final conclusions are presented in Sections 7 
and 8, respectively. 


2 Hybrid methods and order of truncation errors 


Hybrid methods used for solving stiff differential equations are mostly effi- 
cient and have better response than other numerical methods. Consider the 
following initial value problem (IVP): 


eof), eeR”, wig) wa, i <t< ty, (1) 


where f : [to,t¢] x R” — R”. There are several explicit and implicit meth- 
ods for solving such problems, but hybrid methods are more accurate than 
Runge-Kutta and backward differential formulas methods and have a wide 
range of stability; see [7, 5, 6, 4, 10]. Let us consider IVP in the form of (1). 
Linear k-step methods of the following form have 2k +1 arbitrary parameter: 


Cea = Og Peep Oe ee a et Hie + oitece hs 


(2) 
where fait = f (ta + hi tai), fan — f(a — MA, tam) and J, = fla, Tn) 
for m= 1,2,...,k—1. For increasing order of k-step methods in the form of 


(2), a linear combination of the slopes is used at several points between t, 
and ty41, where tn41 = t, +h in which h is the step length on [to, ty]. Then, 
the modified form of (2) with m slopes is given by 


k k m 
Ln+1 = ajay +h>> By fn—j41 +hS" Uj fnsv;s (3) 
j=l j=0 j=l 
where a;, 8;, and yu; are 2k+m-+1 arbitrary parameters; see [10]. Methods 
of form (3) with m off-step points are called hybrid methods, where 
0<v; <1, vj ER, JST Qin Ms 
and herein, we set k = 1 and m = 4. Hence, we write (3) as 


Intl =A1Ly, (4) 
+ h{Bofn+i + Bifa} + h{ i fnter + bafntos + Mafntus + Hafntes}; 
where Q1,Q@2, 30, 61, and v;s are arbitrary parameters and v; is not equal to 


0 or 1 for i = 1, 2, 3, 4. Expanding each term in (4), in Taylor’s series about 
tn, we can obtain a family of the 6th-order methods if the equations 
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Table 1: Values of v;s in cases No. 


1, 2, 3, 4, and 5 of new proposed methods 
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Vi case 1 case 2 case 3 case 4 case 5 
V1 2.020e+00 2.020e+00 2.020e+00 2.020e+00 2.020e+00 
v2 4.000e-01 2.000e-01 4.000e-01 4.000e-01 1.000e-01 
v3 -2.000e-01 -3.000e-01 -2.000e-01 -2.000e-01 -2.000e-01 
va | -1.840e+00 | -1.540e+00 | -1.440e+00 | -1.140e+00 | -1.840e+00 


Table 2: Values of vjs in cases No. 6, 7, 8, 9, and 10 of new proposed methods 
Vi case 6 case 7 case 8 case 9 case 10 
v1 | 2.020e+00 2.090e+00 2.020e+00 | 2.020e+00 | 2.020e+00 
v2 4.000e-01 4.000e-01 3.000e-01 3.000e-01 4.000e-01 
v3 | -2.000e-01 -5.000e-01 -3.000e-01 | -2.000e-01 | -2.0000e-01 
va | -1.840e+00 | -1.840e+00 | -1.840e+00 | -9.540e-02 | -5.240e-01 
ay = 1, 


Bot+ Pr ter +eg+e3+e4 = 1, 

Bo + 101 + Cov2q + ¢303 + C4U4 = 5, 
cyvt + c2v3 + c3v3 + cau) = &, 
a (Bo + cruz + cov} + c3v3 + caf) = 53, 
24 (Bo + evi + e203 + e305 + cavi) Te 
790 (0 + cup + c2v3 + c3¥3 + C4vg) = za; 
759 (Go + crvf + cov$ + c3v§ + cau?) = soap, 


are satisfied, where the principal term of the truncation error is 
1 
ork a™ (tn) + o(h8), cp =1—7(Go + crv; + cous + czvi + cavi). 


In this paper, the coefficients of (4) are proposed and obtained by searching 
and using coefficients that are more stable than the explicit Runge-Kutta 
methods. Values of v; are determined with a wide stability region by searching 
for the interval of [—2.5, 2.5]. The v; values are presented in ten cases as shown 
in Tables 1 and 2. They are called as the new proposed methods in this 
paper: 


Intl = ta Hh Bidet isn, + be fnto. +13 fr+u3 +b tn+os + Bafa fs (5) 


where fn41 = f(tn th, tn41), fato, = f(tnt+vih, fnto,), and fr = f(tn, In) 
for 1 = 1,2,3,4. Note that %41, Un+»;, and x, are numerical approximations 
according to the exact values of the solution x(t) at tr41 = tn +h, trey; = 
ty, + vjh. For converting (5) into explicit methods at each step, the values of 
Ln41 and Lp4,, are predicted and used on right-hand side of the proposed 
methods using the 6th-order explicit Runge-Kutta (RK6) method, respec- 
tively, as follows: 
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Ln+1 = Ln + TRq h(Oky + 64k3 + 49ks5 + 49ke + 9k7), 
ky = Af (fn, Yn); 


ko = hf (an cans hyn + 7 

kg =hf (an + gh Un (3k + kz), 

ka = hf (tn + 3h, Yn + 53 (8k1 + 2k + 8ks), 
ks =hf (tn + 7(7 — 721), Yn + ksn), 


ksn = 399 (3(3V21 — 7)ki — 8(7 — 21) ka + 48(7 — 21)ks) 
+ 5ha(—3(21 ~ V21)ka), (6) 
ke = hf (tn + 4(7 + V21), yn + qo (—5(231 + BW 21) ky 
—40(7 + W21)k2 + ken) 
ken = —320(/21)k3 + 3(21 + 1217/21)k4 + 392(6 + 7/21)ks), 
ky = hf (an +h, Yn + ggg (15(22 + 721) ky + 120k. 
+40(7/21 — 5)k3 + krn), 
krpn = —63(3-721 — 2)k4 — 14(49 + 9V/21)ks + 70(7 — V/21)keo. 


= h 
In41=ln+ 180 — (9k, + 64k3 + 49ks + 49kg + 9k7), (7) 


= uih 
In+u, = * 730 (9k, + 64k3; + 49k; + A9k6; + 9k7i), — 1, 2, 3, A, 


In+1 = tn th{ Pi frotbifnio, csc POST ackas ta Fang POOL ia ts (8) 


where k3;, ks;, ke;, k7i;, 7 = 1,2,3,4, can be obtained by using the method 6 
in which h is replaced by v;h and 


fa4i = Fltn a Py Rata dy frat; = Ihe ais OTs Brabus)» fra= Tbs a). 


Now suppose that the order of (8) is 6 similar to (6). Thus, the difference 
between exact and numerical solutions would be as follows: 


£(tntm) — 2r4m = Crh’2?) (tn) + O(R'). (9) 
The difference operator associated with the 6th-order (8) can be written as 
t(tn+1) — tng. = Ch’ a (t,) + O(R8), 
where C is the error constant of (8). Therefore, we have the following theo- 
rem. 


Theorem 1. Given that (8) is of order p, then, p is equal to 6. 


Proof. Suppose that i = 1,2,3,4 and that x, is exact. According to (7) and 
(8), one can write 
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4 


tat) — &n41 =h)> bi Lt ters £(tatnge)) — TGs Fata,)| 
i=1 


+ hBolf (tna, 2(tn41)) — f(tn41Fn41)] 
+ ChPx)(t,) + O(hPt?). 


Considering the properties of the IVPs in (1), some values such as ny, and 
m belong to intervals of (Zr4v,,2(tn+v,)) and (Fn41, c(tn41)), respectively. 
Thus, we can write 


T Cnvaps@ One) —T Cees Eat) So. li Poy) BC) — Bea) 


F (trois @(tn41)) — f(tmtis Ent) = 92 (tm4i M41) (@(tn41) — En4s)- 


Therefore, using (9), we have 


4 
elt) — fn41 =h S- Li [FE (tnt > Nada) (Elaiee) = Tntv.)| 
i=1 
+ hBo [$e (tna, TMn+1)(2(tn41) — Tn+1)| 
+ Cha) (t,) + O(hPt?). 
Applying (9) to this, we have 


4 
E(tn41) — Un41 =e Mi [FE (tn+osstmntus)Cr ha (tn) + o(ns*)| 


i=l 
+ hBo [$2 (tnt, m+1)Orh a (tn) 7 (ne )] 
+ ChPx) (t,) + O(hP+) 


4 
=H > bi [FE (tntosstmus)Or pre.) | ! 
i=1 


+h® { 8 [32 (tna, tngr)Cyn® +1200) (t,)| 


+Cx!(ty)} + O(R8). 


Thus, it can be concluded that the order of (8) is 6. 


3 Stability analysis of the new methods 


In this section, stability analysis is done on new methods. Dahlquist test 
problem is considered to investigate stability region of the methods presented 
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in this study. Applying the Dahlquist test problem to (5) and inserting p = 6, 
the following equations can be obtained: 


rae 7, (mh)? | (mh)® | (mh)* | (mh)? | (mh)? 
Enum = (Lb mh ry + 3i + Z + 5I a 6! 


(10) 
m= 1, U1, U2, U3, V4, 


where h = hA and pt; € R. Now, let us consider the new hybrid method 
presented in this work: 


In41 = Inth{Bifnt Bo f ntitia f nto, +H2 f nto. tus f n+oz tba f ntvs}- 
(11) 
Substituting (10) into (11), the following equation is obtained: 


_ —  (h)? ey Wale 2 ie 
In+1 tn +h Basin + Bat (1 n+ 4 i 4) 7 5! | 6! )} 
_ ‘ - uh)? jh)? uih)4 
ifm (: r(oh) +! S 7 - [ 
(vjh)® (uih)® 
By yy 


Inserting x, =r" into (11) and dividing it by r”, we can obtain 


prtl — pn {1+ ayh + agh? + agh® ayh? ash’ +apn* 4 azh"}. 


>r=1+aht agh? + agh? + agh* + ash® + agh® + azh’, 


4 4 
a1 =f + fot > mi, a2 = Bot D> (vin), 
i=1 


i=1 


4 4 
a3 = 5 (60+ S-(02m)), a4 = (G0 + S08), 
4= 1. 4=1 


4 


4 
to = py(d0 + So(etH)), a6 = sap (Bo + Doha), 


A 4 
az = 759 (Po + d_ (vi) 


which is a stability polynomial of (11). Figure 1 compares the stability 
region of methods 1 and 2 with that of the 6th-order Runge-Kutta (RK6) 
method (note that method i means that the new method related to the case 
i, 2 = 1,2,3,4). As can be clearly seen, the stability region of the proposed 
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Figure 1: (a) Stability region of the method 1 and RK6 method. (b) Stability region 
of the method 2 and RK6 method. 


jo} 
Oo 


(a) (b) 


Figure 2: (a) Stability region of the method 3 and RK6 method. (b) Stability region 
of the method 4 and RK6 method. 


methods 1 and 2 is much wider than that of the RK6 method. Stability region 
of the proposed methods 3 and 4 is presented in Figure 2 and is compared with 
that of the RK6 method. According to Figure 2, methods 3 and 4 presented 
in this paper have a wider range of stability than the RK6 method, showing 
the efficiency of the proposed methods. Methods 5 and 6 are compared with 
RK6 method in terms of stability region in Figure 3. As demonstrated in 
Figure 3, the proposed methods 5 and 6 have a wider stability region than 
the RK6 method, showing that the efficiency of the proposed methods is 
higher than the RK6 method. 

Figure 4 compares the stability region of methods 7 and 8 with that of the 
RK6 method. As can be seen, methods 7 and 8 have a wider range of stability 
than the RK6 method. Figure 5 also shows the superiority of the methods 
proposed in this paper over the RK6 method. As depicted in Figure 5, the 
width of stability region of methods 9 and 10 is higher compared to the RK6 
method. 
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-- case. 
-~RKG 


~ case.6 
~ PKS 


(b) 
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Figure 3: (a) Stability region of the method 5 and RK6 method. (b) Stability region 


of the method 6 and RK6 method. 


— case7 
“YT | = RKB 


ron 


Figure 4: (a) Stability region of the method 7 and RK6 method. (b) Stability region 


of the method 8 and RK6 method. 


(b) 


Figure 5: (a) Stability region of the method 9 and RK6 method. (b) Stability region 


of the method 10 and RK6 method. 
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4 Optimal control problems 


An OCP includes a cost function J(x,u), a set of state variables, x € X, 
and a set of control variables, u € U. An OCP is solved to find a piecewise 
continuous control u(t), t9 < t < ty, and the associated continuous state 
variable x(t), in order to minimize the given objective function. For more 
explanation, we need a few definitions. 


Definition 1 (Lagrange and Bolza problems). The basic problem in La- 
grange form is 


J(a,u) = / "Gt wa) alO)ak (12) 


to 


Adding another term to functional (12), the Bolza problem is obtained: 
ty 
J(u) = h(tyelt,)) + f° g(t,n(t),ult)at 
to 
An OCP is : 
f 
max J = | g(t, x(t), u(t))dt, (13) 
U to 
a'(t) = f(t,x(d), u(t), 
x(to) =X. 
Note that min{.J} = — max{—/}; see [19]. 


Definition 2 (Hamiltonian). Consider the OCP (13). The function H(t, x, u, ) 
is called as the Hamiltonian function and is equal to 


A, xv, U, Ad) = g(t, xz, u) + Af (t, zr, u) 
where is an adjoint variable. 


Theorem 2 (Pontryagin Maximum Principle). Consider the OCP (13). Sup- 
pose that g(t,xz,u) and f(t,2,u) are both continuously differentiable func- 
tions in their three arguments and concave in « and u. If u* is a control with 
associated state x* and A is a piecewise differentiable function such that u*, 
x*, and X together are satisfied 


A 
ee ee ae 
Ou 
OH 
a Paes 


Ate) = 0, 
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on to <t < ty, then 
J(a*,u*) > J(a,u), 


for any admissible pair (x, w). 


Proof. We refer readers to [13]. 


Theorem 3. [13]. Let the set of controls for problem (13), be Lebesgue 
integrable functions and let tp) < t < ty in R. Suppose that f(t,z,u) is 
concave in u and there exist constants c,,c2,c3 > 0,c4 and @ > 1, such that 


f(t,x,u) = a(t, x) + B(t,x)u, 
|f(@, x, u)| < (1 + |x| + |ul), 
|f(t, 21, u) — f(t,v,u)| < calv1 — 2|(1 + |u)), 
g(t, x,u) < c3|ul® — ca, 


for all t with to) <t<ty , %1,72,uUER. 
Then there exists an optimal pair (x*,u*) maximizing J, with finite J(a*, u*). 


5 Forward-backward sweep method 


For numerically solving the OCP (13) using the indirect method, an algorithm 
is introduced according to the literature [13]. For solving such problems 
numerically first, an algorithm that generates an approximation to an optimal 
piecewise continuous control u*, must divide the time interval of [to, t,] into 
pieces with specific points of interest to = b1,b2,...,bn,6n41 = ty; and 
these points will usually be equally spaced. Approximation will be a vector 
Ud = (u1,U2,---,UN+1), Where u; © u(b;). Any solution to the above OCP 
must also be satisfied: 


nN = 3 A(t) = 0, 
a =0 at uw 


The third equation, namely, optimality condition, can usually be manipu- 
lated to find a representation of u* in terms of t,z, and A. Then, the first 
two equations form a TPBVP. The generalized problem can be solved us- 
ing indirect methods, which are numerical techniques used for solving. The 
FBSM is one of these methods. A rough outline of the algorithm is given 
below. 
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—> 
Here, @ = (@1,%2,.-.,Un41) amd A = (Aq, Ag,...,AN41) are the vector 
approximations for the state and adjoint. 


1. Make an initial guess for UW over the interval. 

2. Using the initial condition 2; = (to) = a and the value for 7, solve @& 
forward in time according to its differential equation in the optimality 
system. 

3. Using the transversality condition Ay+1 = A(t¢) = 0 and the values 

—> 
for @ and @, solve X backward in time according to its differential 


equation in optimality system. 


—> 
4. Update 74 by entering the new @ and X values into the characteriza- 
tion of the optimal control. 


5. Check convergence. If values of variables in this iteration and the last 
iteration are negligibly close, then the current values are considered as 
output solutions. If values are not close, then return to Step 2. 


6 Convergence of FBSM 


For notational simplicity, we express the problem as finding (a(t), A(t), u(t)) 
such that 


a'(t) = f(t, x(t), u(t)), 2(to) = ao, 
N (t) = ki(t, x(t), u(t) + A(t)ka(t, 2(t), u(t), A(t) = 0, 
u(t) = k(t, x(t), u(t). 


Here, to € R” and to < ty are the given real numbers. For a convergence 
analysis of the FBSM, we will make the following assumptions: 
(T) The functions of f,k1, ke, and k3 are Lipschitz continuous with re- 
spect to their second and third arguments, with Lipschitz constants of 
Ly, Ley, Leg, Lez. Moreover, A= |All and H = ||kallo0 < oo. 


Theorem 4. Under the assumptions (T), if 
co =Lg {lexp(L y(t — to)) — Uf + Leg (Le, + ALz2) 
jqlexw(F(ty ~ to) ~ Allexp(y(ty ~to)) +} <1 
then the FBSM is convergent, that is, as n > co, 


— 2”) — y(n) =) 
Oe ee Oa eee 


Proof. We refer the reader to [16]. 
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7 Numerical results 


In this section, examples of various types of OCPs are solved using the pro- 
posed methods, and their numerical results are compared with those of the 
FBSM using the RK method of order 6 (FBSM—RK6). One of the most 
important OCPs is the linear regulator problem, which is generally defined 
as follows. 


Example 1. Let £, Q(t), and R(t) be symmetric and nonnegative definite 
matrices of appropriate dimensions. The so-called linear regulator problem 
(with a linear state-space description) involves a cost functional of the form 


a= 527 (ts) Baty) : ‘ aT Q(t)a(t) + ut (t)R(t)u(t)]|dt. 


For example, we consider an OCP as follows [16]: 


min 5 i (x(t)? + u(t) ]de 


U 


s.t. v'(t) = —a(t) + u(t), x(0) =1. 


The Pontryagin’s maximum principle can be used to construct an analytic 
solution 


H(t,2,u,) = 5 (elt)? + u(t)?) + A(—a(t) + u(t), 


HT 
Ou 
OH 
! 
=~ =- 1)=0. 
N=-Z-=-2+), All)=0 


Together with the state equation, the result will be the following linear dif- 
ferential algebraic system: 


x (t)=—a2(t)+ u(t), 2«(0) =1, 
d (t) =A(t) — a(t), AC) =0, u(t) = —A(2). 


The solution is 


(t) _ V2 cosh(V/2(t — 1)) — sinh(/2(t — 1) 
a V2cosh(/2) + sinh(V2) 
u(t) <—__SimelBUt=) 

V2 cosh(V2) + sinh(V2) 


Optimal value of the objective functional is J = 0.1929092981. The final 
values of the state and optimal are z(1) = 0.2819695346 and 0, respectively, 
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r 
a FBSM-Rk6 
exact il 


() 01 #02 O08 04 O05 O06 O07 O08 09 1 
taxis 


== == FBSM-Rk6 | | 
exact 


u-axes 


0 01 #02 08 04 O05 O06 O07 O08 09 1 
taxis 


(a) 


=—sst= new method-case1 
exact 


x-axis 
° 
fo>) 


oO 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 
taxis 


= = = new method-case1 
exact ‘| 


u-axes 
° 
iy 


10) 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 
taxis 


(b) 


Figure 6: (a) Optimal state and control values of Example 1 FBSM—RK6. (b) Optimal 
state and control values of Example 1 (new proposed method) 


and the initial value of the co-state is (0) = 0.3858185962. Numerical results 
1 


of the problem are shown in Figure 6 with h = 75 and in Tables 3 and 4. 
Numerical results presented in Tables 3 and 4 indicate that each of new ten 
suggested methods calculates the amount of control variable values much 
more accurately than the FBSM—RK6 method. Figure 6 also indicates that 
the new suggested method is exactly based on figure of analytical answer. 
For avoiding overstatement in this paper, one of the diagrams was selected 
and drawn. The rest of figures are similar to each other. Approximate 
performance index is calculated for the proposed method, and the results 
are presented in Tables 5 and 6. According to the results, the precision 
of performance index of new proposed methods is more than that of the 
FBSM—RK6 method by two digits. The Pontryagin’s Theorem is also used 
to solve linear-quadratic problems. 
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Table 3: Error of control values in Example 1 for FBSM—RK6 and new proposed 


Methods 

t h FBSM_RK6 case 1 case 2 case 3 case 4 case 5 
0.90 5 3.0520e-4 7.3292e-5 | 2.5423e-5 | 7.2638e-5 | 7.2172e-5 | 1.4836e-5 
0.90 = 5.3457e-5 2.1141e-5 | 1.5067e-5 | 2.1006e-5 | 2.0916e-5 | 4.1655e-6 
0.90 790 2.6278e-5 1.0962e-5 | 1.4606e-6 | 1.0895e-5 | 1.0850e-5 | 2.5154e-6 
0.90 | sa 1.3036e-5 5.5876e-7 | 8.4673e-7 | 5.5537e-6 | 5.5314e-6 | 1.3740e-6 
Table 4: Error of control values in Example 1 for FBSM—RK6 and new proposed 
methods 

t h FBSM_RK6 case 6 case 7 case 8 case 9 case 10 
0.90 iG 3.0520e-4 6.7372e-5 | 1.1476e-4 | 3.9201le-5 | 2.9350e-5 | 7.1123e-5 
0.90 = 5.3457e-5 1.9984e-5 | 1.5321e-5 | 5.6249e-7 | 1.3363e-6 | 2.0701e-5 
0.90 90 2.6278e-5 1.0386e-5 | 7.1974e-6 | 1.6043e-7 | 1.1054e-6 | 1.0742e-5 
0.90 | a0 1.3036e-5 5.2998e-6 | 3.4747e-6 | 1.9879e-7 | 6.7022e-7 | 5.4774e-6 

Table 5: Errors of the performance index approximation in Examplel 

h FBSM_RK6 case 1 case 2 case 3 case 4 case 5 

a 4.0877e-4 2.8316e-5 | 4.7743e-5 | 2.7735e-5 | 2.7293e-5 | 4.3302e-5 
790 1.9310e-4 6.8092e-6 | 3.0969e-5 | 6.5235e-6 | 6.3012e-6 | 2.8836e-5 
390 9.2656e-5 4.877e-7 1.8338e-5 | 3.4612e-7 | 2.3451le-7 | 1.7294e-5 
7000 1.8014e-5 2.6262e-7 | 4.0178e-6 | 2.9075e-7 | 3.1314e-7 | 3.8123e-6 

Table 6: Errors of the performance index approximation in Example 1 

h FBSM_RK6 case 6 case 7 case 8 case 9 case 10 

5 4.0877e-4 2.3141e-5 | 1.0728e-4 | 6.1400e-5 | 5.4438e-5 | 2.6399e-5 
790 1.9310e-4 4.2415e-6 | 6.0399e-5 | 3.7814e-5 | 3.4359e-5 | 5.8643e-6 
390 9.2656e-5 7.9113e-7 | 3.2969e-5 | 2.1765e-5 | 2.0043e-5 | 1.8678e-8 
a 1.8014e-5 5.1760e-7 | 6.9304e-6 | 4.7036e-6 | 4.3604e-6 | 3.5590e-7 
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Table 7: Control values errors in Example 2 by using FBSM—RK6 and new proposed 


methods 
t h FBSM—RK6 case 1 case 2 case 3 case 4 case 5 
0.70 is 8.2317e-3 1.3560e-4 | 3.7467e-4 | 1.3073e-4 | 1.3148e-4 | 3.1883e-4 
0.70 =a 8.1160e-4 3.7614e-5 | 6.5252e-5 | 3.6716e-5 | 3.6772e-5 | 5.6347e-5 
0.70 T90 4.0487e-4 1.9678e-5 | 3.179le-5 | 1.9236e-5 | 1.9258e-5 | 2.7488e-5 
0.70 | x50 1.0203e-4 1.0221e-5 | 1.5522e-5 | 1.0001le-5 | 1.0010e-5 1.3409-5 


Example 2. Consider the following OCP [18]: 


For solving the above example, using the FBSM and proposed methods, we 


should apply the Pontryagin’s Theorem as follows: 


GG cad) =2n(t)? + 5a(tu(t) | u(t)? | Natt + u(t)), 
aH ; ot 
Ha. =0 at wu >u = —X = 9” 
OH 10 1 1 
ri _ = 
N= ee oe 5s (1) = 0. 


Analytical solutions are as follows [18]: 


ey (tanh(1 — ¢) + .5) cosh(1 — ¢) 
oo cosh(1) : 
«72, _cosh(1 —#) 
as cosh(1) 


The state variable at the end point is (1) = 6.4805427366388¢ — 1. Variable 
control endpoint is u(1) = —3.24027136831le — 1. Optimal value of the ob- 
jective function is J* = 0.3807970779. The proposed methods of Example 2 


were determined as follows in MATLAB environment: 


Table 8: Control values errors in Example 2 by using FBSM—RK6 and new proposed 


methods 
t h FBSM—RK6 case 6 case 7 case 8 case 9 case 10 
0.70 a 8.2317e-3 9.8758e-5 | 8.176le-4 | 4.4379e-4 | 3.9543e-4 | 1.2035e-4 
0.70 sa 8.1160e-4 3.0223e-5 | 5.6542e-5 | 8.1127e-5 | 8.5984e-5 | 3.4750e-5 
0.70 190 4.0487e-4 1.5985e-5 | 7.3752e-5 | 3.9860e-5 | 3.5044e-5 | 1.8260e-5 
0.70 a0 2.0203e-4 8.3745e-6 | 3.6430e-5 | 1.9590e-5 | 1.7183e-5 | 9.5155e-6 
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Figure 7: (a) Optimal state and control values of Example 2 by using FBSM—RK6. 
(b) Optimal state and control values of Example 2 by using new proposed method 


300 Ebadi, Malih Maleki and Ebadian 


Table 9: Errors of the performance index approximation in Example 2 
FBSM—RK6 case 1 case 2 case 3 case 4 case 5 
5.5920e-4 7.6989e-5 | 2.1442e-5 | 7.6139e-5 | 7.6159e-5 | 1.2794e-5 
90 2.9187e-4 2.9588e-5 | 1.9712e-5 | 2.9168e-5 | 2.9173e-5 | 1.5522e-5 
290 1.5038e-4 1.1235e-5 | 1.3437e-5 | 1.1026e-5 | 1.1027e-5 | 1.1376e-5 
3.0605e-5 1.8661le-6 | 3.0719e-6 | 1.8246e-6 | 1.8246e-6 | 2.6652e-6 


HZ]4 


dl 


Numerical results presented in Tables 7 and 8 indicate that each of the new 
ten suggested methods calculates the amount of control variable values much 
more accurately than the FBSM—RK6 method. Figure 7 also indicates that 
the new proposed method is exactly based on figure of analytical answer. For 
simplicity of reporting the results, one of diagrams was selected and drawn. 
The rest of figures are similar to each other. Numerical results presented in 
Tables 9 and 10 indicate that the estimated performance index of the new 
methods is more precise than those of the FBSM—RK6 methods. 


Example 3. Consider the following OCP for a fixed T [14]: 


T ib 
min [Cf e(m)an-+ (u(t))*)at 


s.t. a’ (t)=—a(t)+u(t), 2x(0) =a. 


For converting the problem into the standard form, we can add another state 
and obtain two-dimensional system as follows: 


a(t) = x(t), 


a(t) = i (n)dn, 


We redefine x(t) := [x1(t), v2(t)]7. Thus, we have an OCP of the form: 


Table 10: Errors of the performance index approximation in Example 2 


h FBSM—RK6 case 6 case 7 case 8 case 9 case 10 
= 5.5920e-4 6.9979e-5 | 1.027e5-4 3.6537e-5 2.7287e-5 | 7.4275e-5 
790 2.9187e-4 2.6081e-5 | 6.0153e-5 2.7391e-5 2.2764e-5 | 2.8242e-5 
290 1.5038e-4 9.4814e-6 | 3.3603e-5 | 1.73105e-5 | 1.4995e-5 | 1.0565e-5 
Tog 3.0605e-5 1.5153e-6 | 7.0964e-6 3.8519e-6 3.3889e-6 | 1.7326e-6 
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Figure 8: (a) Optimal state and control values of Example 3 using FBSM—RK6. (b) 
Optimal state and control values of Example 3 using new proposed method 


T 
min | (xo (t) + u(t)2)dt 


s.t., x(t) =—a(t) + u(t), 


Analytical solution of the problem is 


A =(x2 + u’) + Ai(—21 + u) + A241, 


A 1 
° =2u+ A; = 0 at ut = ul = —5A1, A(L) = A2(T) = 0, 
U 
OH OH 
ee ey. Me ee 
> 1 Ox, 1 2) Ax> 4 
=\() =—(t -T) —1 + et), 
u*(t) = — S(t) = (1 +4-T- et), 


Numerical results for Example 3 are obtained and shown in Figure 8 and 
Tables 11 and 12. 


Table 11: Control values errors in Example 3 using FBSM—RK6 and new proposed 


method 
t h FBSM—RK6 case 1 case 2 case 3 case 4 case 5 
0.80 5 1.3027e-3 4.3609e-5 | 1.9320e-4 | 4.5033e-5 | 4.5546e-5 | 1.7794e-4 
0.80 =5 2.4376e-4 1.6351e-6 | 2.9366e-5 | 1.8584e-6 | 1.9768e-6 | 2.7439e-5 
0.80 T90 1.2086e-4 4.1195e-7 | 1.4147e-5 | 2.2016e-7 | 5.8037e-7 | 1.3244e-5 
0.80 | a55 6.0174e-5 1.0924e-7 | 6.9448e-6 | 1.625le-7 | 1.9287e-7 | 6.5116e-6 
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Table 12: Control values errors in Example 3 by using FBSM—RK6 and new proposed 


methods 
t h FBSM—RK6 case 6 case 7 case 8 case 9 case 10 
0.80 5 1.3027e-3 5.4085e-5 | 3.2411e-4 | 2.1473e-4 | 2.0072e-4 | 4.8041e-5 
0.80 a 2.4376e-4 3.5454e-6 | 5.1626e-5 | 3.4084e-5 | 3.1524e-5 | 2.3610e-6 
0.80 T90 1.2086e-4 1.3566e-6 | 2.5050e-5 | 1.6529e-5 | 1.5263e-5 | 7.6572e-7 
0.80 | 50 6.0174e-5 5.7899e-7 | 1.2340e-5 | 8.1415e-6 | 7.5122e-6 | 2.8398¢-7 


Numerical results presented in Tables 11 and 12 indicate that each of the 
new ten suggested methods calculates the amount of control variable values 
much more accurately than the FBSM—RK6 method. Figure 8 also indicates 
that the figure of the new proposed method is quite matched on real answer 
and is much better than the FBSM—RK6 method. For simplicity of reporting 
the results, one of diagrams was selected and drawn. The rest of figures are 
similar to each other. 


8 Conclusion 


A new class of the 6th-order explicit hybrid methods was presented for which 
the 6th-order Runge-Kutta method is used as a predictor scheme to gain 
whole method of the same order, and the order of truncation errors was in- 
vestigated for the explicit hybrid Runge-Kutta methods. The stability of 
the methods was discussed, and the results revealed that the stability re- 
gions of the proposed methods are wider compared to the 6th-order explicit 
Runge-Kutta method. Finally, three examples of OCPs were solved using 
MATLAB, FBSM scheme, and the presented methods and numerical results 
related to given examples were presented in Tables 3-12. According to the 
findings, it can be concluded that the new explicit hybrid methods have a 
good performance in accuracy and performance index approximation com- 
pared to the RK6 method. 
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